﻿using System;
using System.Collections.Generic;
using System.Globalization;
using System.Linq;
using System.Text;

namespace IMPression
{
	class ParseException : Exception
	{
		public string ParsedString {get; set;}
		public int IndexOfError {get; set;}
		public ParseException(string parsedString, int indexOfError, string message) : base(message)
		{
			ParsedString = parsedString;
			IndexOfError = indexOfError;
		}
	}
	public abstract class EquationElement
	{
		public EquationElement Parent {get; set;}
		protected StringBuilder m_sb = new StringBuilder();
		public virtual int Length { get { return m_sb.Length; }}
		protected virtual void Add(char ch)
		{
			m_sb.Append(ch);
		}
		public abstract bool IsValid(string equation, int index);
		public override string ToString()
		{
			return m_sb.ToString();
		}
		protected EquationElement Root
		{
			get
			{
				if (Parent == null)
					return this;
				var parent = Parent;
				while (parent != null && parent.Parent != null)
					parent = parent.Parent;
				return parent;
			}
		}

		Dictionary<string, double> m_variables;
		Dictionary<string, bool> m_usedVars;
		public void ClearVars()
		{
			var root = Root;
			if (root.m_variables != null)
				root.m_variables.Clear();
		}
		public bool VarExist(string varname)
		{
			var root = Root;
			return (root.m_variables != null && root.m_variables.ContainsKey(varname.ToLower()));
		}
		public void SetVar(string varname, string value)
		{
			var term = new Term();
			term.Parse(value);
			SetVar(varname, term.Value);
		}
		public void SetVar(string varname, double value)
		{
			var root = Root;
			if (root.m_variables == null)
				root.m_variables = new Dictionary<string,double>();
			root.m_variables[varname.ToLower()] = value;
		}
		public double GetVar(string varname)
		{
			var root = Root;
			double result = 0;
			if (root.m_variables == null || root.m_variables.TryGetValue(varname.ToLower(), out result) == false)
				throw new ParseException("", 0, string.Format("La variable '{0}' n'est pas définie", varname));
			if (m_usedVars == null)
				m_usedVars = new Dictionary<string, bool>();
			if (m_usedVars.ContainsKey(varname) == false)
				m_usedVars[varname] = true;
			return result;
		}
		public string[] GetUsedVars()
		{
			if (m_usedVars == null)
				return new string[0];
			var l = new List<string>(m_usedVars.Keys);
			l.Sort();
			return l.ToArray();
		}
	}
	class Operator : EquationElement
	{
		static public bool IsValidOperator(char ch)
		{
			var operators = new[] { '+', '-', '/', '*', '^', '%' };
			return operators.Contains(ch);
		}
		public char Value { get; private set;}
		public override bool IsValid(string equation, int index)
		{
			return IsValidOperator(equation[index]);
		}
		public Operator(string src, int start)
		{
			if (IsValid(src, start))
			{
				Value = src[start];
				Add(Value);
			}
		}
		public Operator(char ch)
		{
			if (IsValidOperator(ch))
			{
				Value = ch;
				Add(Value);
			}
		}
		public bool CanStartTerm
		{
			get { return Value == '-'; }
		}
		public bool IsSign
		{
			get { return Value == '-'; }
		}
		public double Perform(double left, double right)
		{
			if (Value == '-')
				return left - right;
			if (Value == '+')
				return left + right;
			if (Value == '*')
				return left * right;
			if (Value == '/')
				return left / right;
			if (Value == '^')
				return Function.Pow(left, right);
			if (Value == '%')
				return left % right;
			throw new ParseException("", 0, string.Format("Unsupported operator '{0}'", Value));
		}
	}
	public abstract class EquationValue : EquationElement
	{
		protected bool m_signed;
		public virtual double Value {get; protected set;}
		public virtual bool Signed 
		{
			get { return m_signed;}
			set
			{
				m_signed = value;
				if (m_signed && Value > 0)
				{
					Value = 0 - Value;
					return;
				}
				if (!m_signed && Value < 0)
				{
					Value = Math.Abs(Value);
				}

			}
		}
		public virtual bool CanBeSigned { get {return true;} }
		public override string ToString()
		{
			if (Signed)
				return "s" + base.ToString();
			return base.ToString();
		}
	}
	class Term : EquationValue
	{
		public override double Value
		{
			get
			{
				double value =	0;
				if (Stack.Count == 0)
					return value;

				var left = Stack[0] as EquationValue;
				value = left.Value;

				for (var index = 1; index < Stack.Count; index++)
				{
					Operator op = null;
					EquationValue right = null;
					if (index < Stack.Count)
						op = Stack[index] as Operator;
					if (index+1 < Stack.Count)
						right = Stack[index+1] as EquationValue;
					if (op != null && right != null)
						value = op.Perform(value, right.Value);
					index += 1;
				}
				if (Signed)
					return -value;
				return value;
			}
		}
		public List<EquationElement> Stack { get; set; }
		static public bool IsValidDigit(char ch)
		{
			if (ch == '(' || ch == '[')
				return true;
			return false;
		}
		public override bool IsValid(string equation, int index)
		{
			return IsValidDigit(equation[index]);
		}
		public Term()
		{
			Stack = new List<EquationElement>();
		}
		public override string ToString()
		{
			var sb = new StringBuilder();
			if (Signed)
				sb.Append("s");
			sb.Append("(");
			foreach (var item in Stack)
				sb.Append(item);
			sb.Append(")");
			return sb.ToString();
		}
		public void Parse(string equation)
		{
			Parse(equation, this);
		}
		public void Parse(string equation, EquationElement root)
		{
			Stack.Clear();
			Parse(equation, 0, root);
			CombineTerms(new[] {'^'});
			CombineTopDown();
			CombineTerms(new[] {'*', '/', '%'});
		}
		EquationElement PopSignOperator()
		{
			if (Stack.Count > 1)
			{
				var index = Stack.Count-1;
				var op1 = Stack[index] as Operator;
				var op2 = Stack[index-1] as Operator;
				if (op2 != null && op1 != null && op1.IsSign)
				{
					Stack.RemoveAt(index);
					return op1;
				}
			}
			if (Stack.Count == 1)
			{
				var index = Stack.Count-1;
				var op1 = Stack[index] as Operator;
				if (op1 != null && op1.IsSign)
				{
					Stack.RemoveAt(index);
					return op1;
				}
			}
			return null;
		}
		EquationElement LastElement
		{
			get
			{
				if (Stack.Count == 0)
					return null;
				return Stack[Stack.Count-1];
			}
		}
		int OperatorCount
		{
			get
			{
				var cnt = 0;
				for (var index = Stack.Count-1; index >= 0; index--)
				{
					if (Stack[index] is Operator)
					{
						cnt++;
						continue;
					}
					break;
				}
				return cnt;
			}
		}
		EquationElement Add(EquationElement item)
		{
			item.Parent = this;
			if (item is EquationValue)
			{
				// insérer l'opérateur '*' entre les valeurs
				if (LastElement is EquationValue) 
					Stack.Add(new Operator('*'));
				// si il y a 2 opérateurs et que le dernier est -, le supprimer et marquer la valeur comme signée (négative)
				if (PopSignOperator() != null)
					((EquationValue)item).Signed = true;
				Stack.Add(item);
				return item;
			}
			if (item is Operator)
			{
				// seul - peut démarrer un terme
				if (((Operator)item).CanStartTerm == false && Stack.Count == 0)
					throw new ParseException(string.Empty, 0, "Term cannot start with operator");
				// si il y a un second opérateur s'assurer que c'est un signe
				var op = LastElement as Operator;
				if (op != null && ((Operator)item).IsSign == false)
					throw new ParseException("", 0, "consecutive operators are not valid");
				if (OperatorCount >= 2)
					throw new ParseException("", 0, "consecutive operators are not valid");
				Stack.Add(item);
				return item;
			}
			return null;
		}
		int Parse(string equation, int index, EquationElement root)
		{
			while (index < equation.Length)
			{
				var ch = equation[index];
				if (char.IsWhiteSpace(ch))
				{
					index++;
					continue;
				}
				if (Operator.IsValidOperator(ch))
				{
					var n = Add(new Operator(equation, index));
					index += n.Length;
					continue;
				}
				if (Number.IsValidDigit(equation, index))
				{
					var n = Add(new Number(equation, index));
					index += n.Length;
					continue;
				}
				if (Function.IsValidDigit(equation, index))
				{
					var n = Add(new Function(equation, index, root));
					index += n.Length;
					continue;
				}
				if (Constant.IsValidDigit(equation, index))
				{
					var n = Add(new Constant(equation, index));
					index += n.Length;
					continue;
				}				
				if (Variable.IsValidDigit(equation, index, root))
				{
					var n = Add(new Variable(equation, index, root));
					index += n.Length;
					continue;
				}
				index++;

				if (IsValidDigit(ch))
				{
					var endindex = FindMatchingEnd(ch, equation, index-1);
					if (endindex > index)
					{
						var len = endindex - index;
						var s = equation.Substring(index, len);
						var g = Add(new Term()) as Term;
						len = g.Parse(s, 0, root) + 1;
						index += len;
						continue;
					}
					throw new ParseException(equation, index-1, "Pas de parenthèses fermantes correspondantes trouvées");
				}

			}
			return index;
		}
		void CombineTopDown()
		{
			// traitement spécial de l'opérateur de puissance. Si un opérateur est élevé à la puissance d'une puissance, la puissance doit être calculée en premier.
			// ex. 2^2^3 = 2^(2^3) = 256. Ceci est réalisé en combinant de droite à gauche (ou de haut en bas)
			var operatorCh = '^';
			foreach (var element in Stack)
			{
				if (element is Term)
					((Term)element).CombineTopDown();
			}
			var index = Stack.Count-1;
			while (index > 0)
			{
				var op = Stack[index] as Operator;
				index--;
				if (op == null || op.Value != operatorCh)
					continue;

				var left = Stack[index];
				var right = Stack[index+2];

				var newterm = new Term {Parent = this};
				newterm.Add(left);
				newterm.Add(op);
				newterm.Add(right);

				if ( ((EquationValue)left).Signed )
				{
					((EquationValue)left).Signed = false;
					newterm.Signed = true;
				}

				Stack.RemoveRange(index, newterm.Stack.Count);
				Stack.Insert(index, newterm);
			}
		}
		void CombineTerms(char[] operators)
		{
			foreach (var element in Stack)
			{
				if (element is Term)
					((Term)element).CombineTerms(operators);
			}
			if (NeedSubTerms(operators) == false)
				return;

			var startIndex = 0;
			while (startIndex < Stack.Count)
			{
				startIndex = FindOperator(startIndex, operators);
				if (startIndex < 0)
					return;

				var newterm = new Term {Parent = this};
				startIndex--;
				var startpoint = startIndex;

				while (startIndex < Stack.Count)
				{
					var item = Stack[startIndex];
					if (item is EquationValue)
					{
						newterm.Add(item); 
						startIndex++;
					}
					if (item is Operator)
					{
						var op = item as Operator;
						if (op == null || operators.Contains(op.Value) == false)
						{
							Stack.RemoveRange(startpoint, newterm.Stack.Count);
							Stack.Insert(startpoint, newterm);
							break;
						}
						newterm.Add(item);
						startIndex++;
					}

					if (startIndex >= Stack.Count)
					{
						Stack.RemoveRange(startpoint, newterm.Stack.Count);
						Stack.Insert(startpoint, newterm);
						return;
					}
				}
			}
		}
		bool NeedSubTerms(IEnumerable<char> operators)
		{
			// si tous les opérateurs d'un terme sont du type correct, alors pas besoin de sous-termes.
			return Stack.OfType<Operator>().Any(op => operators.Contains(op.Value) == false);
		}

		int FindOperator(int startIndex, char[] operators)
		{
			for (var index = startIndex; index < Stack.Count; index++)
			{
				var op = Stack[index] as Operator;
				if (op != null && operators.Contains(op.Value))
					return index;
			}
			return -1;
		}

		public static int FindMatchingEnd(char beginChar, string equation, int beginCharindex)
		{
			var index = beginCharindex;
			var matchCount = 0;
			var endChar = ')';
			if (beginChar == '[')
				endChar = ']';
			while (index < equation.Length-1)
			{
				index++;
				var ch = equation[index];
				if (ch == beginChar)
				{
					matchCount++;
					continue;
				}
				if (ch == endChar)
				{
					if (matchCount == 0)
						return index;
					matchCount--;
				}
			}
			return -1;
		}
		public static string ExtractName(string equation, int index)
		{
			return ExtractName(equation, index, true);
		}
		public static string ExtractName(string equation, int index, bool allowDigits)
		{
			var sb = new StringBuilder();
			while (index < equation.Length)
			{
				var ch = equation[index];
				if (char.IsLetter(ch))
					sb.Append(ch);
				else if (allowDigits && (char.IsDigit(ch)) && (sb.Length > 0))  // chiffre autorisé si il n'est pas en première position
					sb.Append(ch);
				else
					break;
				index++;
			}
			return sb.ToString();
		}
	}
	class Number : EquationValue
	{
		static string GetNumber(string equation, int index)
		{
			var sb = new StringBuilder();
			var chars = new[] {',', '.', 'e', 'E','-','+', '0', '1', '2', '3', '4', '5', '6', '7', '8', '9'};
			var isScientific = false;
			while (index < equation.Length)
			{
				var ch = equation[index++];
				if (ch >= '0' && ch <= '9')
				{
					sb.Append(ch);
					continue;
				}
				if (chars.Contains(ch))
				{
					// support pour l'écriture scientifique
					if (ch == 'e' || ch == 'E')
					{
						if (isScientific)
							break;
						isScientific = true;
						if (index >= equation.Length)
							break;
						var nextchar = equation[index];
						if (chars.Contains(nextchar) == false)
							break;
					}
					// +- seulement autorisé si précédé par 'e'
					var lastchar = char.ToLower(sb[sb.Length-1]);
					if (ch == '-' && lastchar != 'e')
						break;
					if (ch == '+' && lastchar != 'e')
						break;
					sb.Append(ch);
					continue;
				}
				break;
			}
			return sb.ToString();
		}
		static public bool IsValidDigit(string equation, int index)
		{
			var ch = equation[index];
			if (ch < '0' || ch > '9')
				return false;
			return GetNumber(equation, index).Length > 0;
		}
		public override bool IsValid(string equation, int index)
		{
			return IsValidDigit(equation, index);
		}
		public Number(string src, int start)
		{
			var number = GetNumber(src, start);
			if (number.Length > 0)
				m_sb.Append(number);
			double result;

			var val = m_sb.ToString();
			var invarSeperator = CultureInfo.InvariantCulture.NumberFormat.NumberDecimalSeparator;
			val = val.Replace(",", invarSeperator);
			if (double.TryParse(val, NumberStyles.Float, CultureInfo.InvariantCulture, out result) == false)
				throw new ParseException(src, start, "Impossible de convertir en nombre à virgule");
			Value = result;
		}
	}   
	class Variable : EquationValue
	{
		public override bool Signed 
		{
			get { return m_signed;}
			set { m_signed = value; }
		}
		public override double Value
		{
			get
			{
				var value = Root.GetVar(m_sb.ToString());
				if (Signed)
					return -value;
				return value;
			}
		}
		static string GetVarName(string equation, int index, EquationElement root)
		{
			var varname = Term.ExtractName(equation, index).ToLower();
			if (root != null && root.VarExist(varname))
				return varname;
			var ch = equation[index];
			if (ch >= 'a' && ch <= 'z')
				return new string(new[] {ch});
			if (ch >= 'A' && ch <= 'Z')
				return new string(new[] {ch});
			return string.Empty;
		}
		static public bool IsValidDigit(string equation, int index, EquationElement root)
		{
			return GetVarName(equation, index, root).Length > 0;
		}
		public override bool IsValid(string equation, int index)
		{
			throw new NotImplementedException();
		}
		public bool IsValid(string equation, int index, EquationElement root)
		{
			return GetVarName(equation, index, root).Length > 0;
		}
		public Variable(string src, int start, EquationElement root)
		{
			// variable d'un seul caractère supportée pour le moment :/
			var varname = GetVarName(src, start, root);
			if (varname.Length > 0)
				m_sb.Append(varname);
		}
	}
	public class Constant : EquationValue
	{
		static bool IsPI(string equation, int index)
		{
			if (equation.Length > index + 1)
			{
				var spi = Term.ExtractName(equation, index, false).ToLower();
				if (spi == "pi")
					return true;
			}
			return false;
		}
		static bool IsPhi(string equation, int index)
		{
			if (equation.Length > index + 1)
			{
				var spi = Term.ExtractName(equation, index, false).ToLower();
				if (spi == "phi")
					return true;
			}
			return false;
		}
		static bool IsEulerGamma(string equation, int index)
		{
			if (equation.Length > index + 1)
			{
				var spi = Term.ExtractName(equation, index, false).ToLower();
				if (spi == "eulergamma")
					return true;
			}
			return false;
		}
        static bool IsCatalan(string equation, int index)
        {
            if (equation.Length > index + 1)
            {
                var spi = Term.ExtractName(equation, index, false).ToLower();
                if (spi == "catalan")
                    return true;
            }
            return false;
        }
		static bool IsE(string equation, int index)
		{
			var ch = equation[index];
			if (ch == 'e' || ch == 'E')
				return true;
			return false;
		}
        static bool IsSpeedOfLight(string equation, int index)
        {
            var ch = equation[index];
            if (ch == 'c' || ch == 'C')
                return true;
            return false;
        }
		static public bool IsValidDigit(string equation, int index)
		{
		    return IsE(equation, index) || IsPI(equation, index) || IsPhi(equation, index) || IsEulerGamma(equation, index) ||
		           IsCatalan(equation, index) || IsSpeedOfLight(equation, index);
		}

	    public override bool IsValid(string equation, int index)
		{
			return IsValidDigit(equation, index);
		}
        public static double E { get { return 2.7182818284590452353602874713526624977572470937000d; } }

        public static double Pi { get { return 3.1415926535897932384626433832795028841971693993751d; } }

        public static double Phi { get { return 1.6180339887498948482045868343656381177203091798057628621354486227052604628189024497072d; } }

        public static double EulerGamma { get { return 0.5772156649015328606065120900824024310421593359399235988057672348849d; } }

        public static double Catalan { get { return 0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694794d; } }

        public static double SpeedOfLight { get { return 2.99792458e8; } }

		public Constant(string src, int start)
		{
			if (IsE(src, start))
			{
				m_sb.Append("E");
				Value = E;
                return;
			}
			if (IsPI(src, start))
			{
				m_sb.Append("pi");
				Value = Pi;
                return;
			}
			if(IsPhi(src, start))
			{
				m_sb.Append("phi");
				Value = Phi;
                return;
			}
			if(IsEulerGamma(src, start))
			{
				m_sb.Append("eulergamma");
				Value = EulerGamma;
                return;
			}
            if(IsCatalan(src, start))
            {
                m_sb.Append("catalan");
                Value = Catalan;
                return;
            }
            if(IsSpeedOfLight(src, start))
            {
                m_sb.Append("c");
                Value = SpeedOfLight;
            }
		}
	}
	public class Function : EquationValue
	{
		delegate double FunctionCallbackNone();
		delegate double FunctionCallback(double value);
		delegate double FunctionCallback2(double value, double value1);
        delegate double FunctionCallback3(double value, double value1, double value2);
        delegate double FunctionCallback4(double value, double value1, double value2, double value3);
        delegate double FunctionCallbackInf(params double[] values);
		public static List<string> FunctionsList { get; private set; }
		static Dictionary<string, FunctionCallback> Functions;
		static Dictionary<string, FunctionCallback2> Functions2;
        static Dictionary<string, FunctionCallback3> Functions3;
        static Dictionary<string, FunctionCallback4> Functions4;
        static Dictionary<string, FunctionCallbackInf> FunctionsInf;
		static Dictionary<string, FunctionCallbackNone> FunctionsNone;

		private static Random _rnd = new Random();

	    private static bool degrees
	    {
	        get { return EquationParser.UseDegrees; }
	    }

	    private static double DegToRad(double deg)
        {
            return (Constant.Pi / 180) * deg;
        }
        private static double RadToDeg(double rad)
        {
            return rad * (180.0 / Constant.Pi);
        }

        #region Fonctions trigonométriques
        public static double Sec(double a)
        {
            return 1 / Math.Cos(degrees ? DegToRad(a) : a);
        }
        public static double Csc(double a)
        {
            return 1 / Math.Sin(degrees ? DegToRad(a) : a);
        }
        public static double Cot(double a)
        {
            return 1 / Math.Tan(degrees ? DegToRad(a) : a);
        }
        public static double Sech(double a)
        {
            return 1 / Math.Cosh(degrees ? DegToRad(a) : a);
        }
        public static double Csch(double a)
        {
            return 1 / Math.Sinh(degrees ? DegToRad(a) : a);
        }
        public static double Coth(double a)
        {
            return 1 / Math.Tanh(degrees ? DegToRad(a) : a);
        }
        public static double Arcsec(double a)
        {
            return (degrees ? RadToDeg(Math.Acos(1 / a)) : Math.Acos((1 / a)));
        }
		public static double Arccsc(double a)
		{
            return (degrees ? RadToDeg(Math.Acos(1 / a)) : Math.Asin((1 / a)));
		}
		public static double Arccot(double a)
		{
            return (degrees ? RadToDeg(Constant.Pi / 2 - Math.Atan(a)) : Constant.Pi / 2 - Math.Atan(a));
		}

		public static double Arcsinh(double a)
		{
		    return (degrees ? RadToDeg(Ln(a + Math.Sqrt(a * a + 1))) : Ln(a + Math.Sqrt(a * a + 1)));
		}
		public static double Arccosh(double a)
		{
		    return (degrees ? RadToDeg(Ln(a + Math.Sqrt(a * Math.Sqrt(a - 1)))) : Ln(a + Math.Sqrt(a * Math.Sqrt(a - 1))));
		}
		public static double Arctanh(double a)
		{
		    return (degrees
		        ? RadToDeg(Ln(1 / a + Math.Sqrt(1 / a + 1) * Math.Sqrt(1 / a - 1)))
		        : Ln(1 / a + Math.Sqrt(1 / a + 1) * Math.Sqrt(1 / a - 1)));;
		}
		public static double Arccoth(double a)
		{
            return (degrees ? RadToDeg(0.5 * Ln((a + 1) / (a - 1))) : 0.5 * Ln((a + 1) / (a - 1)));
		}
		public static double Arccsch(double a)
		{
            return (degrees ? RadToDeg(Ln(1 / a + Math.Sqrt(1 / (a * a + 1)))) : Ln(1 / a + Math.Sqrt(1 / (a * a + 1))));
		}
		public static double Arcsech(double a)
		{
            return (degrees ? RadToDeg(Ln(1 / a + Math.Sqrt(1 / a + 1) * Math.Sqrt(1 / a - 1))) : Ln(1 / a + Math.Sqrt(1 / a + 1) * Math.Sqrt(1 / a - 1)));
		}

		public static double Sinc(double a)
		{
            return Math.Sin(degrees ? DegToRad(a) : a) / (degrees ? DegToRad(a) : a);
		}

        public static double Cos(double a)
        {
            return Math.Cos(degrees ? DegToRad(a) : a);
        }
        public static double Sin(double a)
        {
            return Math.Sin(degrees ? DegToRad(a) : a);
        }
        public static double Tan(double a)
        {
            return Math.Tan(degrees ? DegToRad(a) : a);
        }
        public static double Cosh(double a)
        {
            return Math.Cosh(degrees ? DegToRad(a) : a);
        }
        public static double Sinh(double a)
        {
            return Math.Sinh(degrees ? DegToRad(a) : a);
        }
        public static double Tanh(double a)
        {
            return Math.Tanh(degrees ? DegToRad(a) : a);
        }
        #endregion

        #region Fonctions entières
        public static double Ln(double value)
        {
            return Math.Log(value, Constant.E);
        }
        public static double Rand()
        {
            return _rnd.NextDouble();
        }
        public static double RandN(double a, double b)
        {
            if (b <= a)
            {
                double c = a;
                a = b;
                b = c;
            }
            return _rnd.Next((int)a, (int)b);
        }
        public static double Rand(double a, double b)
        {
            if (b <= a)
            {
                double c = a;
                a = b;
                b = c;
            }
            return _rnd.NextDouble() * (b - a) + a;
        }
        public static double Pow(double x, double d)
        {
            if (Math.Truncate(d) == d)
            {
                var n = (int) (Math.Truncate(d));
                if (n < 0) return (1.0 / Math.Pow(x, -n));
                switch (n)
                {
                    case 0:
                        return (1.0);
                    case 1:
                        return (x);
                    case 2:
                        return (x * x);
                    case 3:
                        return (x * x * x);
                    case 4:
                    {
                        double x2 = x * x;
                        return (x2 * x2);
                    }
                    case 5:
                    {
                        double x2 = x * x;
                        return (x2 * x2 * x);
                    }
                    case 6:
                    {
                        double x2 = x * x;
                        return (x2 * x2 * x2);
                    }
                    case 7:
                    {
                        double x3 = x * x * x;
                        return (x3 * x3 * x);
                    }
                    case 8:
                    {
                        double x2 = x * x;
                        double x4 = x2 * x2;
                        return (x4 * x4);
                    }
                    case 9:
                    {
                        double x3 = x * x * x;
                        return (x3 * x3 * x3);
                    }
                    case 10:
                    {
                        double x2 = x * x;
                        double x4 = x2 * x2;
                        return (x4 * x4 * x2);
                    }
                    case 12:
                    {
                        double x3 = x * x * x;
                        double x6 = x3 * x3;
                        return (x6 * x6);
                    }
                    case 16:
                    {
                        double x2 = x * x;
                        double x4 = x2 * x2;
                        double x8 = x4 * x4;
                        return (x8 * x8);
                    }
                    default:
                        return (Math.Pow(x, n));
                }
            }
            else
            {
                return Math.Pow(x, d);
            }
        }
        public static double Eq(double a, double b)
        {
            return a == b ? 1 : 0;
        }
        public static double Fact(double a)
        {
            if (Math.Truncate(a) != a || a < 0) return double.NaN;
            return a > 1 ? a * Fact(a - 1) : 1;
        }
        public static double GCD(double a, double b)
        {
            return b == 0 ? a : GCD(b, a % b);
        }
        public static double GCD(params double[] a)
        {
            return a.Aggregate(GCD);
        }
        public static double LCM(double a, double b)
        {
            return (a * b) / GCD(a, b);
        }
        public static double LCM(params double[] a)
        {
            return a.Aggregate(LCM);
        }
        public static double Fibo(double a)
        {
            //return (1 / Math.Sqrt(5)) * (Math.Pow(Constant.Phi, a) - Math.Pow(-(1 / Constant.Phi), a));
            return (Pow(Constant.Phi, a) - Math.Cos(a * Constant.Pi) * Pow(Constant.Phi, -a)) / Math.Sqrt(5);
        }
        public static double Fact2(double a)
        {
            return Math.Pow(2 / Constant.Pi, (0.25 * (1 - Math.Cos(Constant.Pi * a)))) * Math.Pow(2, (a / 2)) * Gamma(a / 2 + 1);
        }
        public static double Round(double a, double b)
        {
            return Math.Round(a, (int)b);
        }
        public static double CubeRoot(double a)
        {
            return (Math.Pow(a, (double)1 / 3));
        }
        public static double NthRoot(double a, double n)
        {
            return Math.Pow(a, (double)1 / n);
        }
        public static double Max(params double[] a)
        {
            return a.Max();
        }
        public static double Min(params double[] a)
        {
            return a.Min();
        }

        public static double KroneckerDelta(params double[] a)
        {
            return a.All(x => x == a[0]) ? 1 : 0;
        }
        public static double KroneckerDelta(double a)
        {
            return a == 0 ? 1 : 0;
        }
        public static double Average(params double[] a)
        {
            return a.Sum() / (double)a.Count();
        }
        public static double LucasL(double a)
        {
            return Pow(Constant.Phi, a) + Pow(Constant.Phi, -a) * Math.Cos(a * Constant.Pi);
        }
        public static double ProductLog(double a)
        {
            return ProductLog(0, a);
        }
        public static double ProductLog(double a, double z)
        {
            double runningSum = 0;
            for (int n = 1; n <= (int)z; ++n)
            {
                runningSum += Pow(-n, n - 1) / Fact(n) * Pow(a, n);
            }
            return runningSum;
        }
        public static double Mediane(params double[] a)
        {
            if (a.Length == 0) return double.PositiveInfinity;
            if (a.Length == 1) return a[0];
            Array.Sort(a);
            if(a.Length % 2 == 0)
            {
                return Average(a[a.Length / 2 - 1], a[a.Length / 2]);
            }
            else
            {
                return a[(a.Length - 1) / 2];
            }
        }
        public static double Q1(params double[] a)
        {
            if (a.Length == 0) return double.PositiveInfinity;
            if (a.Length == 1) return a[0];
            Array.Sort(a);
            return a[(int)Math.Ceiling((double)(a.Length - 1) / 4.0) - 1];
        }
        public static double Q3(params double[] a)
        {
            if (a.Length == 0) return double.PositiveInfinity;
            if (a.Length == 1) return a[0];
            Array.Sort(a);
            return a[(int)Math.Ceiling((double)(a.Length - 1) * (3.0/4.0)) - 1];
        }
        public static double Somme(params double[] a)
        {
            return a.Sum();
        }
        #endregion

        #region Fonctions Beta
        public static double BetaLn(double z, double w)
        {
            if (z <= 0.0)
            {
                throw new ParseException(z.ToString(), 0, "z doit être positif.");
            }
            if (w <= 0.0)
            {
                throw new ParseException(w.ToString(), 0, "w doit être positif.");
            }
            return GammaLn(z) + GammaLn(w) - GammaLn(z + w);
        }
        public static double Beta(double z, double w)
        {
            return Math.Exp(BetaLn(z, w));
        }
        public static double BetaIncomplete(double a, double b, double x)
        {
            return BetaRegularized(a, b, x) * Beta(a, b);
        }
        public static double BetaRegularized(double a, double b, double x)
        {
            if (a < 0.0)
            {
                throw new ParseException(a.ToString(), 0, "a doit être positif.");
            }
            if (b < 0.0)
            {
                throw new ParseException(b.ToString(), 0, "b doit être positif.");
            }
            if (x < 0.0 || x > 1.0)
            {
                throw new ParseException(x.ToString(), 0, "x doit être entre 0 et 1.");
            }
            var bt = (x == 0.0 || x == 1.0)
                ? 0.0
                : Math.Exp(GammaLn(a + b) - GammaLn(a) - GammaLn(b) + (a * Math.Log(x)) + (b * Math.Log(1.0 - x)));
            var symmetryTransformation = x >= (a + 1.0) / (a + b + 2.0);
            var eps = Math.Pow(2, -53);
            var fpmin = 1.0 / eps;
            if (symmetryTransformation)
            {
                x = 1.0 - x;
                var swap = a;
                a = b;
                b = swap;
            }
            var qab = a + b;
            var qap = a + 1.0;
            var qam = a - 1.0;
            var c = 1.0;
            var d = 1.0 - (qab * x / qap);
            if (Math.Abs(d) < fpmin)
            {
                d = fpmin;
            }
            d = 1.0 / d;
            var h = d;
            for (int m = 1, m2 = 2; m <= 140; m++, m2 += 2)
            {
                var aa = m * (b - m) * x / ((qam + m2) * (a + m2));
                d = 1.0 + (aa * d);
                if (Math.Abs(d) < fpmin)
                {
                    d = fpmin;
                }
                c = 1.0 + (aa / c);
                if (Math.Abs(c) < fpmin)
                {
                    c = fpmin;
                }
                d = 1.0 / d;
                h *= d * c;
                aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
                d = 1.0 + (aa * d);
                if (Math.Abs(d) < fpmin)
                {
                    d = fpmin;
                }
                c = 1.0 + (aa / c);
                if (Math.Abs(c) < fpmin)
                {
                    c = fpmin;
                }
                d = 1.0 / d;
                var del = d * c;
                h *= del;
                if (Math.Abs(del - 1.0) <= eps)
                {
                    return symmetryTransformation ? 1.0 - (bt * h / a) : bt * h / a;
                }
            }
            return symmetryTransformation ? 1.0 - (bt * h / a) : bt * h / a;
        }
        #endregion

        #region Fonctions Gamma
        const int GammaN = 10;
        const double GammaR = 10.900511;
        public static double Gamma(double x)
        {
            if (x <= 0.0)
            {
                string msg = string.Format("Argument invalide {0}. L'argument doit être positif.", x);
                throw new ParseException(x.ToString(), 0, msg);
            }
            const double gamma = 0.577215664901532860606512090;
            if (x < 0.001)
                return 1.0 / (x * (1.0 + gamma * x));
            if (x < 12.0)
            {
                double y = x;
                int n = 0;
                bool arg_was_less_than_one = (y < 1.0);
                if (arg_was_less_than_one)
                {
                    y += 1.0;
                }
                else
                {
                    n = (int)(Math.Floor(y)) - 1;
                    y -= n;
                }
                double[] p =
				{
					-1.71618513886549492533811E+0,
					 2.47656508055759199108314E+1,
					-3.79804256470945635097577E+2,
					 6.29331155312818442661052E+2,
					 8.66966202790413211295064E+2,
					-3.14512729688483675254357E+4,
					-3.61444134186911729807069E+4,
					 6.64561438202405440627855E+4
				};
                double[] q =
				{
					-3.08402300119738975254353E+1,
					 3.15350626979604161529144E+2,
					-1.01515636749021914166146E+3,
					-3.10777167157231109440444E+3,
					 2.25381184209801510330112E+4,
					 4.75584627752788110767815E+3,
					-1.34659959864969306392456E+5,
					-1.15132259675553483497211E+5
				};
                double num = 0.0;
                double den = 1.0;
                int i;
                double z = y - 1;
                for (i = 0; i < 8; i++)
                {
                    num = (num + p[i]) * z;
                    den = den * z + q[i];
                }
                double result = num / den + 1.0;
                if (arg_was_less_than_one)
                {
                    result /= (y - 1.0);
                }
                else
                {
                    for (i = 0; i < n; i++)
                        result *= y++;
                }
                return result;
            }
            if (x > 171.624)
            {
                return double.PositiveInfinity;
            }
            return Math.Exp(LogGamma(x));
        }

        public static double PolyGamma(double z)
        {
            double x = 0;
            double k = 1;

            while (true)
            {
                double nw = 1 / k - 1 / (k + z - 1);
                if(Math.Round(x + nw, 15) == Math.Round(x, 15))
                {
                    break;
                }
                x += nw;

                k++;
            }

            return x - Constant.EulerGamma;
        }

        public static double LogGamma(double x)
        {
            if (x <= 0.0)
            {
                string msg = string.Format("Argument invalide {0}. L'argument doit être positif.", x);
                throw new ParseException(x.ToString(), 0, msg);
            }
            if (x < 12.0)
            {
                return Math.Log(Math.Abs(Gamma(x)));
            }
            double[] c =
			{
				 1.0/12.0,
				-1.0/360.0,
				1.0/1260.0,
				-1.0/1680.0,
				1.0/1188.0,
				-691.0/360360.0,
				1.0/156.0,
				-3617.0/122400.0
			};
            double z = 1.0 / (x * x);
            double sum = c[7];
            for (int i = 6; i >= 0; i--)
            {
                sum *= z;
                sum += c[i];
            }
            double series = sum / x;
            double halfLogTwoPi = 0.91893853320467274178032973640562;
            double logGamma = (x - 0.5) * Math.Log(x) - x + halfLogTwoPi + series;
            return logGamma;
        }
        static readonly double[] GammaDk =
        {
            2.48574089138753565546e-5,
            1.05142378581721974210,
            -3.45687097222016235469,
            4.51227709466894823700,
            -2.98285225323576655721,
            1.05639711577126713077,
            -1.95428773191645869583e-1,
            1.70970543404441224307e-2,
            -5.71926117404305781283e-4,
            4.63399473359905636708e-6,
            -2.71994908488607703910e-9
        };
        public static double GammaLn(double z)
        {
            if (z < 0.5)
            {
                double s = GammaDk[0];
                for (int i = 1; i <= GammaN; i++)
                {
                    s += GammaDk[i] / (i - z);
                }

                return 1.1447298858494001741434273513530587116472948129153d
                       - Math.Log(Math.Sin(Math.PI * z))
                       - Math.Log(s)
                       - 0.6207822376352452223455184457816472122518527279025978
                       - ((0.5 - z) * Math.Log((0.5 - z + GammaR) / Math.E));
            }
            else
            {
                double s = GammaDk[0];
                for (int i = 1; i <= GammaN; i++)
                {
                    s += GammaDk[i] / (z + i - 1.0);
                }

                return Math.Log(s)
                       + 0.6207822376352452223455184457816472122518527279025978
                       + ((z - 0.5) * Math.Log((z - 0.5 + GammaR) / Math.E));
            }
        }
        // buggé
        /*public static double Gamma(double z)
        {
            if (z < 0.5)
            {
                double s = GammaDk[0];
                for (int i = 1; i <= GammaN; i++)
                {
                    s += GammaDk[i] / (i - z);
                }
                return Math.PI / (Math.Sin(Math.PI * z)
                                * s
                                * 1.8603827342052657173362492472666631120594218414085755
                                * Math.Pow((0.5 - z + GammaR) / Math.E, 0.5 - z));
            }
            else
            {
                double s = GammaDk[0];
                for (int i = 1; i <= GammaN; i++)
                {
                    s += GammaDk[i] / (z + i - 1.0);
                }

                return s * 1.8603827342052657173362492472666631120594218414085755 * Math.Pow((z - 0.5 + GammaR) / Math.E, z - 0.5);
            }
        }*/
        public static double GammaUpperRegularized(double a, double x)
        {
            const double epsilon = 0.000000000000001;
            const double big = 4503599627370496.0;
            const double bigInv = 2.22044604925031308085e-16;
            if (x <= 0d || a <= 0d)
            {
                return 1d;
            }
            if (x < 1d || x < a)
            {
                return 1d - GammaLowerRegularized(a, x);
            }
            double ax = a * Math.Log(x) - x - GammaLn(a);
            if (ax < -709.78271289338399)
            {
                return a < x ? 0d : 1d;
            }
            ax = Math.Exp(ax);
            double t;
            double y = 1 - a;
            double z = x + y + 1;
            double c = 0;
            double pkm2 = 1;
            double qkm2 = x;
            double pkm1 = x + 1;
            double qkm1 = z * x;
            double ans = pkm1 / qkm1;
            do
            {
                c = c + 1;
                y = y + 1;
                z = z + 2;
                double yc = y * c;
                double pk = pkm1 * z - pkm2 * yc;
                double qk = qkm1 * z - qkm2 * yc;
                if (qk != 0)
                {
                    double r = pk / qk;
                    t = Math.Abs((ans - r) / r);
                    ans = r;
                }
                else
                {
                    t = 1;
                }
                pkm2 = pkm1;
                pkm1 = pk;
                qkm2 = qkm1;
                qkm1 = qk;
                if (Math.Abs(pk) > big)
                {
                    pkm2 = pkm2 * bigInv;
                    pkm1 = pkm1 * bigInv;
                    qkm2 = qkm2 * bigInv;
                    qkm1 = qkm1 * bigInv;
                }
            }
            while (t > epsilon);
            return ans * ax;
        }
        public static double GammaUpperIncomplete(double a, double x)
        {
            return GammaUpperRegularized(a, x) * Gamma(a);
        }
        public static double GammaLowerIncomplete(double a, double x)
        {
            return GammaLowerRegularized(a, x) * Gamma(a);
        }
        public static double GammaLowerRegularized(double a, double x)
        {
            const double epsilon = 0.000000000000001;
            const double big = 4503599627370496.0;
            const double bigInv = 2.22044604925031308085e-16;
            if (a < 0d)
            {
                throw new ParseException(a.ToString(), 0, "a doit être positif.");
            }
            if (x < 0d)
            {
                throw new ParseException(x.ToString(), 0, "x doit être positif.");
            }
            if (a == 0.0)
            {
                if (x == 0.0)
                {
                    return double.NaN;
                }
                return 1d;
            }
            if (x == 0.0)
            {
                return 0d;
            }
            double ax = (a * Math.Log(x)) - x - GammaLn(a);
            if (ax < -709.78271289338399)
            {
                return a < x ? 1d : 0d;
            }
            if (x <= 1 || x <= a)
            {
                double r2 = a;
                double c2 = 1;
                double ans2 = 1;
                do
                {
                    r2 = r2 + 1;
                    c2 = c2 * x / r2;
                    ans2 += c2;
                }
                while ((c2 / ans2) > epsilon);
                return Math.Exp(ax) * ans2 / a;
            }
            int c = 0;
            double y = 1 - a;
            double z = x + y + 1;
            double p3 = 1;
            double q3 = x;
            double p2 = x + 1;
            double q2 = z * x;
            double ans = p2 / q2;
            double error;
            do
            {
                c++;
                y += 1;
                z += 2;
                double yc = y * c;
                double p = (p2 * z) - (p3 * yc);
                double q = (q2 * z) - (q3 * yc);
                if (q != 0)
                {
                    double nextans = p / q;
                    error = Math.Abs((ans - nextans) / nextans);
                    ans = nextans;
                }
                else
                {
                    error = 1;
                }
                p3 = p2;
                p2 = p;
                q3 = q2;
                q2 = q;
                if (Math.Abs(p) > big)
                {
                    p3 *= bigInv;
                    p2 *= bigInv;
                    q3 *= bigInv;
                    q2 *= bigInv;
                }
            }
            while (error > epsilon);
            return 1d - (Math.Exp(ax) * ans);
        }
        public static double GammaLowerRegularizedInv(double a, double y0)
        {
            const double epsilon = 0.000000000000001;
            const double big = 4503599627370496.0;
            const double threshold = 5 * epsilon;
            if (double.IsNaN(a) || double.IsNaN(y0))
            {
                return double.NaN;
            }
            if (a < 0 || a == 0.0)
            {
                throw new ArgumentOutOfRangeException("a");
            }
            if (y0 < 0 || y0 > 1)
            {
                throw new ArgumentOutOfRangeException("y0");
            }
            if (y0 == 1.0)
            {
                return 0d;
            }
            if (y0 == 1.0)
            {
                return double.PositiveInfinity;
            }
            y0 = 1 - y0;
            double xUpper = big;
            double xLower = 0;
            double yUpper = 1;
            double yLower = 0;
            double d = 1 / (9 * a);
            double y = 1 - d - (0.98 * 1.4142135623730950488016887242096980785696718753769d * ErfInv((2.0 * y0) - 1.0) * Math.Sqrt(d));
            double x = a * y * y * y;
            double lgm = GammaLn(a);
            for (int i = 0; i < 10; i++)
            {
                if (x < xLower || x > xUpper)
                {
                    d = 0.0625;
                    break;
                }
                y = 1 - GammaLowerRegularized(a, x);
                if (y < yLower || y > yUpper)
                {
                    d = 0.0625;
                    break;
                }
                if (y < y0)
                {
                    xUpper = x;
                    yLower = y;
                }
                else
                {
                    xLower = x;
                    yUpper = y;
                }
                d = ((a - 1) * Math.Log(x)) - x - lgm;
                if (d < -709.78271289338399)
                {
                    d = 0.0625;
                    break;
                }
                d = -Math.Exp(d);
                d = (y - y0) / d;
                if (Math.Abs(d / x) < epsilon)
                {
                    return x;
                }
                if ((d > (x / 4)) && (y0 < 0.05))
                {
                    d = x / 10;
                }
                x -= d;
            }
            if (xUpper == big)
            {
                if (x <= 0)
                {
                    x = 1;
                }
                while (xUpper == big)
                {
                    x = (1 + d) * x;
                    y = 1 - GammaLowerRegularized(a, x);
                    if (y < y0)
                    {
                        xUpper = x;
                        yLower = y;
                        break;
                    }
                    d = d + d;
                }
            }
            int dir = 0;
            d = 0.5;
            for (int i = 0; i < 400; i++)
            {
                x = xLower + (d * (xUpper - xLower));
                y = 1 - GammaLowerRegularized(a, x);
                lgm = (xUpper - xLower) / (xLower + xUpper);
                if (Math.Abs(lgm) < threshold)
                {
                    return x;
                }
                lgm = (y - y0) / y0;
                if (Math.Abs(lgm) < threshold)
                {
                    return x;
                }
                if (x <= 0d)
                {
                    return 0d;
                }
                if (y >= y0)
                {
                    xLower = x;
                    yUpper = y;
                    if (dir < 0)
                    {
                        dir = 0;
                        d = 0.5;
                    }
                    else
                    {
                        if (dir > 1)
                        {
                            d = (0.5 * d) + 0.5;
                        }
                        else
                        {
                            d = (y0 - yLower) / (yUpper - yLower);
                        }
                    }
                    dir = dir + 1;
                }
                else
                {
                    xUpper = x;
                    yLower = y;
                    if (dir > 0)
                    {
                        dir = 0;
                        d = 0.5;
                    }
                    else
                    {
                        if (dir < -1)
                        {
                            d = 0.5 * d;
                        }
                        else
                        {
                            d = (y0 - yLower) / (yUpper - yLower);
                        }
                    }
                    dir = dir - 1;
                }
            }
            return x;
        }
        public static double DiGamma(double x)
        {
            const double c = 12.0;
            const double d1 = -0.57721566490153286;
            const double d2 = 1.6449340668482264365;
            const double s = 1e-6;
            const double s3 = 1.0 / 12.0;
            const double s4 = 1.0 / 120.0;
            const double s5 = 1.0 / 252.0;
            const double s6 = 1.0 / 240.0;
            const double s7 = 1.0 / 132.0;
            if (double.IsNegativeInfinity(x) || double.IsNaN(x))
            {
                return double.NaN;
            }
            if (x <= 0 && Math.Floor(x) == x)
            {
                return double.NegativeInfinity;
            }
            if (x < 0)
            {
                return DiGamma(1.0 - x) + (Math.PI / Math.Tan(-Math.PI * x));
            }
            if (x <= s)
            {
                return d1 - (1 / x) + (d2 * x);
            }
            double result = 0;
            while (x < c)
            {
                result -= 1 / x;
                x++;
            }
            if (x >= c)
            {
                var r = 1 / x;
                result += Math.Log(x) - (0.5 * r);
                r *= r;

                result -= r * (s3 - (r * (s4 - (r * (s5 - (r * (s6 - (r * s7))))))));
            }
            return result;
        }
        public static double DiGammaInv(double p)
        {
            if (double.IsNaN(p))
            {
                return double.NaN;
            }
            if (double.IsNegativeInfinity(p))
            {
                return 0.0;
            }
            if (double.IsPositiveInfinity(p))
            {
                return double.PositiveInfinity;
            }
            var x = Math.Exp(p);
            for (var d = 1.0; d > 1.0e-15; d /= 2.0)
            {
                x += d * Math.Sign(p - DiGamma(x));
            }
            return x;
        }
        #endregion

        #region Fonctions Erf
        public static double Erf(double x)
        {
            if (x == 0)
            {
                return 0;
            }
            if (double.IsPositiveInfinity(x))
            {
                return 1;
            }
            if (double.IsNegativeInfinity(x))
            {
                return -1;
            }
            if (double.IsNaN(x))
            {
                return double.NaN;
            }

            return ErfImp(x, false);
        }
        public static double Erfc(double x)
        {
            if (x == 0)
            {
                return 1;
            }

            if (double.IsPositiveInfinity(x))
            {
                return 0;
            }

            if (double.IsNegativeInfinity(x))
            {
                return 2;
            }

            if (double.IsNaN(x))
            {
                return double.NaN;
            }

            return ErfImp(x, true);
        }
        public static double ErfInv(double z)
        {
            if (z == 0.0)
            {
                return 0.0;
            }

            if (z >= 1.0)
            {
                return double.PositiveInfinity;
            }

            if (z <= -1.0)
            {
                return double.NegativeInfinity;
            }

            double p, q, s;
            if (z < 0)
            {
                p = -z;
                q = 1 - p;
                s = -1;
            }
            else
            {
                p = z;
                q = 1 - z;
                s = 1;
            }

            return ErfInvImpl(p, q, s);
        }
        static double ErfImp(double z, double invert)
        {
            if (invert != 1 && invert != 0) return 0;
            return ErfImp(z, invert == 1);
        }
        static double ErfImp(double z, bool invert)
        {
            if (z < 0)
            {
                if (!invert)
                {
                    return -ErfImp(-z, false);
                }

                if (z < -0.5)
                {
                    return 2 - ErfImp(-z, true);
                }

                return 1 + ErfImp(-z, false);
            }

            double result;

            // Big bunch of selection statements now to pick which
            // implementation to use, try to put most likely options
            // first:
            if (z < 0.5)
            {
                // We're going to calculate erf:
                if (z < 1e-10)
                {
                    result = (z * 1.125) + (z * 0.003379167095512573896158903121545171688);
                }
                else
                {
                    // Worst case absolute error found: 6.688618532e-21
                    double[] nc = { 0.00337916709551257388990745, -0.00073695653048167948530905, -0.374732337392919607868241, 0.0817442448733587196071743, -0.0421089319936548595203468, 0.0070165709512095756344528, -0.00495091255982435110337458, 0.000871646599037922480317225 };
                    double[] dc = { 1, -0.218088218087924645390535, 0.412542972725442099083918, -0.0841891147873106755410271, 0.0655338856400241519690695, -0.0120019604454941768171266, 0.00408165558926174048329689, -0.000615900721557769691924509 };

                    result = (z * 1.125) + (z * Polynomial(z, nc) / Polynomial(z, dc));
                }
            }
            else if ((z < 110) || ((z < 110) && invert))
            {
                // We'll be calculating erfc:
                invert = !invert;
                double r, b;
                if (z < 0.75)
                {
                    // Worst case absolute error found: 5.582813374e-21
                    double[] nc = { -0.0361790390718262471360258, 0.292251883444882683221149, 0.281447041797604512774415, 0.125610208862766947294894, 0.0274135028268930549240776, 0.00250839672168065762786937 };
                    double[] dc = { 1, 1.8545005897903486499845, 1.43575803037831418074962, 0.582827658753036572454135, 0.124810476932949746447682, 0.0113724176546353285778481 };
                    r = Polynomial(z - 0.5, nc) / Polynomial(z - 0.5, dc);
                    b = 0.3440242112F;
                }
                else if (z < 1.25)
                {
                    // Worst case absolute error found: 4.01854729e-21
                    double[] nc = { -0.0397876892611136856954425, 0.153165212467878293257683, 0.191260295600936245503129, 0.10276327061989304213645, 0.029637090615738836726027, 0.0046093486780275489468812, 0.000307607820348680180548455 };
                    double[] dc = { 1, 1.95520072987627704987886, 1.64762317199384860109595, 0.768238607022126250082483, 0.209793185936509782784315, 0.0319569316899913392596356, 0.00213363160895785378615014 };
                    r = Polynomial(z - 0.75, nc) / Polynomial(z - 0.75, dc);
                    b = 0.419990927F;
                }
                else if (z < 2.25)
                {
                    // Worst case absolute error found: 2.866005373e-21
                    double[] nc = { -0.0300838560557949717328341, 0.0538578829844454508530552, 0.0726211541651914182692959, 0.0367628469888049348429018, 0.00964629015572527529605267, 0.00133453480075291076745275, 0.778087599782504251917881e-4 };
                    double[] dc = { 1, 1.75967098147167528287343, 1.32883571437961120556307, 0.552528596508757581287907, 0.133793056941332861912279, 0.0179509645176280768640766, 0.00104712440019937356634038, -0.106640381820357337177643e-7 };
                    r = Polynomial(z - 1.25, nc) / Polynomial(z - 1.25, dc);
                    b = 0.4898625016F;
                }
                else if (z < 3.5)
                {
                    // Worst case absolute error found: 1.045355789e-21
                    double[] nc = { -0.0117907570137227847827732, 0.014262132090538809896674, 0.0202234435902960820020765, 0.00930668299990432009042239, 0.00213357802422065994322516, 0.00025022987386460102395382, 0.120534912219588189822126e-4 };
                    double[] dc = { 1, 1.50376225203620482047419, 0.965397786204462896346934, 0.339265230476796681555511, 0.0689740649541569716897427, 0.00771060262491768307365526, 0.000371421101531069302990367 };
                    r = Polynomial(z - 2.25, nc) / Polynomial(z - 2.25, dc);
                    b = 0.5317370892F;
                }
                else if (z < 5.25)
                {
                    // Worst case absolute error found: 8.300028706e-22
                    double[] nc = { -0.00546954795538729307482955, 0.00404190278731707110245394, 0.0054963369553161170521356, 0.00212616472603945399437862, 0.000394984014495083900689956, 0.365565477064442377259271e-4, 0.135485897109932323253786e-5 };
                    double[] dc = { 1, 1.21019697773630784832251, 0.620914668221143886601045, 0.173038430661142762569515, 0.0276550813773432047594539, 0.00240625974424309709745382, 0.891811817251336577241006e-4, -0.465528836283382684461025e-11 };
                    r = Polynomial(z - 3.5, nc) / Polynomial(z - 3.5, dc);
                    b = 0.5489973426F;
                }
                else if (z < 8)
                {
                    // Worst case absolute error found: 1.700157534e-21
                    double[] nc = { -0.00270722535905778347999196, 0.0013187563425029400461378, 0.00119925933261002333923989, 0.00027849619811344664248235, 0.267822988218331849989363e-4, 0.923043672315028197865066e-6 };
                    double[] dc = { 1, 0.814632808543141591118279, 0.268901665856299542168425, 0.0449877216103041118694989, 0.00381759663320248459168994, 0.000131571897888596914350697, 0.404815359675764138445257e-11 };
                    r = Polynomial(z - 5.25, nc) / Polynomial(z - 5.25, dc);
                    b = 0.5571740866F;
                }
                else if (z < 11.5)
                {
                    // Worst case absolute error found: 3.002278011e-22
                    double[] nc = { -0.00109946720691742196814323, 0.000406425442750422675169153, 0.000274499489416900707787024, 0.465293770646659383436343e-4, 0.320955425395767463401993e-5, 0.778286018145020892261936e-7 };
                    double[] dc = { 1, 0.588173710611846046373373, 0.139363331289409746077541, 0.0166329340417083678763028, 0.00100023921310234908642639, 0.24254837521587225125068e-4 };
                    r = Polynomial(z - 8, nc) / Polynomial(z - 8, dc);
                    b = 0.5609807968F;
                }
                else if (z < 17)
                {
                    // Worst case absolute error found: 6.741114695e-21
                    double[] nc = { -0.00056907993601094962855594, 0.000169498540373762264416984, 0.518472354581100890120501e-4, 0.382819312231928859704678e-5, 0.824989931281894431781794e-7 };
                    double[] dc = { 1, 0.339637250051139347430323, 0.043472647870310663055044, 0.00248549335224637114641629, 0.535633305337152900549536e-4, -0.117490944405459578783846e-12 };
                    r = Polynomial(z - 11.5, nc) / Polynomial(z - 11.5, dc);
                    b = 0.5626493692F;
                }
                else if (z < 24)
                {
                    // Worst case absolute error found: 7.802346984e-22
                    double[] nc = { -0.000241313599483991337479091, 0.574224975202501512365975e-4, 0.115998962927383778460557e-4, 0.581762134402593739370875e-6, 0.853971555085673614607418e-8 };
                    double[] dc = { 1, 0.233044138299687841018015, 0.0204186940546440312625597, 0.000797185647564398289151125, 0.117019281670172327758019e-4 };
                    r = Polynomial(z - 17, nc) / Polynomial(z - 17, dc);
                    b = 0.5634598136F;
                }
                else if (z < 38)
                {
                    // Worst case absolute error found: 2.414228989e-22
                    double[] nc = { -0.000146674699277760365803642, 0.162666552112280519955647e-4, 0.269116248509165239294897e-5, 0.979584479468091935086972e-7, 0.101994647625723465722285e-8 };
                    double[] dc = { 1, 0.165907812944847226546036, 0.0103361716191505884359634, 0.000286593026373868366935721, 0.298401570840900340874568e-5 };
                    r = Polynomial(z - 24, nc) / Polynomial(z - 24, dc);
                    b = 0.5638477802F;
                }
                else if (z < 60)
                {
                    // Worst case absolute error found: 5.896543869e-24
                    double[] nc = { -0.583905797629771786720406e-4, 0.412510325105496173512992e-5, 0.431790922420250949096906e-6, 0.993365155590013193345569e-8, 0.653480510020104699270084e-10 };
                    double[] dc = { 1, 0.105077086072039915406159, 0.00414278428675475620830226, 0.726338754644523769144108e-4, 0.477818471047398785369849e-6 };
                    r = Polynomial(z - 38, nc) / Polynomial(z - 38, dc);
                    b = 0.5640528202F;
                }
                else if (z < 85)
                {
                    // Worst case absolute error found: 3.080612264e-21
                    double[] nc = { -0.196457797609229579459841e-4, 0.157243887666800692441195e-5, 0.543902511192700878690335e-7, 0.317472492369117710852685e-9 };
                    double[] dc = { 1, 0.052803989240957632204885, 0.000926876069151753290378112, 0.541011723226630257077328e-5, 0.535093845803642394908747e-15 };
                    r = Polynomial(z - 60, nc) / Polynomial(z - 60, dc);
                    b = 0.5641309023F;
                }
                else
                {
                    // Worst case absolute error found: 8.094633491e-22
                    double[] nc = { -0.789224703978722689089794e-5, 0.622088451660986955124162e-6, 0.145728445676882396797184e-7, 0.603715505542715364529243e-10 };
                    double[] dc = { 1, 0.0375328846356293715248719, 0.000467919535974625308126054, 0.193847039275845656900547e-5 };
                    r = Polynomial(z - 85, nc) / Polynomial(z - 85, dc);
                    b = 0.5641584396F;
                }

                double g = Math.Exp(-z * z) / z;
                result = (g * b) + (g * r);
            }
            else
            {
                // Any value of z larger than 28 will underflow to zero:
                result = 0;
                invert = !invert;
            }

            if (invert)
            {
                result = 1 - result;
            }

            return result;
        }
        static double ErfInvImpl(double p, double q, double s)
        {
            double result;

            if (p <= 0.5)
            {
                // Evaluate inverse erf using the rational approximation:
                //
                // x = p(p+10)(Y+R(p))
                //
                // Where Y is a constant, and R(p) is optimized for a low
                // absolute error compared to |Y|.
                //
                // double: Max error found: 2.001849e-18
                // long double: Max error found: 1.017064e-20
                // Maximum Deviation Found (actual error term at infinite precision) 8.030e-21
                const float y = 0.0891314744949340820313f;
                double[] nc = { -0.000508781949658280665617, -0.00836874819741736770379, 0.0334806625409744615033, -0.0126926147662974029034, -0.0365637971411762664006, 0.0219878681111168899165, 0.00822687874676915743155, -0.00538772965071242932965 };
                double[] dc = { 1, -0.970005043303290640362, -1.56574558234175846809, 1.56221558398423026363, 0.662328840472002992063, -0.71228902341542847553, -0.0527396382340099713954, 0.0795283687341571680018, -0.00233393759374190016776, 0.000886216390456424707504 };
                double g = p * (p + 10);
                double r = Polynomial(p, nc) / Polynomial(p, dc);
                result = (g * y) + (g * r);
            }
            else if (q >= 0.25)
            {
                // Rational approximation for 0.5 > q >= 0.25
                //
                // x = sqrt(-2*log(q)) / (Y + R(q))
                //
                // Where Y is a constant, and R(q) is optimized for a low
                // absolute error compared to Y.
                //
                // double : Max error found: 7.403372e-17
                // long double : Max error found: 6.084616e-20
                // Maximum Deviation Found (error term) 4.811e-20
                const float y = 2.249481201171875f;
                double[] nc = { -0.202433508355938759655, 0.105264680699391713268, 8.37050328343119927838, 17.6447298408374015486, -18.8510648058714251895, -44.6382324441786960818, 17.445385985570866523, 21.1294655448340526258, -3.67192254707729348546 };
                double[] dc = { 1, 6.24264124854247537712, 3.9713437953343869095, -28.6608180499800029974, -20.1432634680485188801, 48.5609213108739935468, 10.8268667355460159008, -22.6436933413139721736, 1.72114765761200282724 };
                double g = Math.Sqrt(-2 * Math.Log(q));
                double xs = q - 0.25;
                double r = Polynomial(xs, nc) / Polynomial(xs, dc);
                result = g / (y + r);
            }
            else
            {
                // For q < 0.25 we have a series of rational approximations all
                // of the general form:
                //
                // let: x = sqrt(-log(q))
                //
                // Then the result is given by:
                //
                // x(Y+R(x-B))
                //
                // where Y is a constant, B is the lowest value of x for which
                // the approximation is valid, and R(x-B) is optimized for a low
                // absolute error compared to Y.
                //
                // Note that almost all code will really go through the first
                // or maybe second approximation.  After than we're dealing with very
                // small input values indeed: 80 and 128 bit long double's go all the
                // way down to ~ 1e-5000 so the "tail" is rather long...
                double x = Math.Sqrt(-Math.Log(q));
                if (x < 3)
                {
                    // Max error found: 1.089051e-20
                    const float y = 0.807220458984375f;
                    double[] nc = { -0.131102781679951906451, -0.163794047193317060787, 0.117030156341995252019, 0.387079738972604337464, 0.337785538912035898924, 0.142869534408157156766, 0.0290157910005329060432, 0.00214558995388805277169, -0.679465575181126350155e-6, 0.285225331782217055858e-7, -0.681149956853776992068e-9 };
                    double[] dc = { 1, 3.46625407242567245975, 5.38168345707006855425, 4.77846592945843778382, 2.59301921623620271374, 0.848854343457902036425, 0.152264338295331783612, 0.01105924229346489121 };
                    double xs = x - 1.125;
                    double r = Polynomial(xs, nc) / Polynomial(xs, dc);
                    result = (y * x) + (r * x);
                }
                else if (x < 6)
                {
                    // Max error found: 8.389174e-21
                    const float y = 0.93995571136474609375f;
                    double[] nc = { -0.0350353787183177984712, -0.00222426529213447927281, 0.0185573306514231072324, 0.00950804701325919603619, 0.00187123492819559223345, 0.000157544617424960554631, 0.460469890584317994083e-5, -0.230404776911882601748e-9, 0.266339227425782031962e-11 };
                    double[] dc = { 1, 1.3653349817554063097, 0.762059164553623404043, 0.220091105764131249824, 0.0341589143670947727934, 0.00263861676657015992959, 0.764675292302794483503e-4 };
                    double xs = x - 3;
                    double r = Polynomial(xs, nc) / Polynomial(xs, dc);
                    result = (y * x) + (r * x);
                }
                else if (x < 18)
                {
                    // Max error found: 1.481312e-19
                    const float y = 0.98362827301025390625f;
                    double[] nc = { -0.0167431005076633737133, -0.00112951438745580278863, 0.00105628862152492910091, 0.000209386317487588078668, 0.149624783758342370182e-4, 0.449696789927706453732e-6, 0.462596163522878599135e-8, -0.281128735628831791805e-13, 0.99055709973310326855e-16 };
                    double[] dc = { 1, 0.591429344886417493481, 0.138151865749083321638, 0.0160746087093676504695, 0.000964011807005165528527, 0.275335474764726041141e-4, 0.282243172016108031869e-6 };
                    double xs = x - 6;
                    double r = Polynomial(xs, nc) / Polynomial(xs, dc);
                    result = (y * x) + (r * x);
                }
                else if (x < 44)
                {
                    // Max error found: 5.697761e-20
                    const float y = 0.99714565277099609375f;
                    double[] nc = { -0.0024978212791898131227, -0.779190719229053954292e-5, 0.254723037413027451751e-4, 0.162397777342510920873e-5, 0.396341011304801168516e-7, 0.411632831190944208473e-9, 0.145596286718675035587e-11, -0.116765012397184275695e-17 };
                    double[] dc = { 1, 0.207123112214422517181, 0.0169410838120975906478, 0.000690538265622684595676, 0.145007359818232637924e-4, 0.144437756628144157666e-6, 0.509761276599778486139e-9 };
                    double xs = x - 18;
                    double r = Polynomial(xs, nc) / Polynomial(xs, dc);
                    result = (y * x) + (r * x);
                }
                else
                {
                    // Max error found: 1.279746e-20
                    const float y = 0.99941349029541015625f;
                    double[] nc = { -0.000539042911019078575891, -0.28398759004727721098e-6, 0.899465114892291446442e-6, 0.229345859265920864296e-7, 0.225561444863500149219e-9, 0.947846627503022684216e-12, 0.135880130108924861008e-14, -0.348890393399948882918e-21 };
                    double[] dc = { 1, 0.0845746234001899436914, 0.00282092984726264681981, 0.468292921940894236786e-4, 0.399968812193862100054e-6, 0.161809290887904476097e-8, 0.231558608310259605225e-11 };
                    double xs = x - 44;
                    double r = Polynomial(xs, nc) / Polynomial(xs, dc);
                    result = (y * x) + (r * x);
                }
            }

            return s * result;
        }
        #endregion

        #region Fonctions polynomiales
        public static double Polynomial(params double[] args)
        {
            double z = args[0];
            double[] coefficients = args;
            var a = coefficients.ToList();
            a.RemoveAt(0);
            coefficients = a.ToArray();
            double sum = coefficients[coefficients.Length - 1];
            for (int i = coefficients.Length - 2; i >= 0; --i)
            {
                sum *= z;
                sum += coefficients[i];
            }

            return sum;
        }

        public static double Polynomial (double z, double[] coeffs)
        {
            var c = coeffs.ToList();
            c.Insert(0, z);
            return Polynomial(c.ToArray());
        }

        public static double HermiteH(double n, double z)
        {
            if(n < 0)
            {
                throw new ParseException("", 0, "n doit être positif.");
            }
            if(n == 0)
            {
                return 1.0;
            }

            double x = 0;

            for (int k = 0; k <= Math.Truncate(n / 2); k++)
            {
                x += (Pow(-1, k) * Pow(2 * z, n - 2 * x)) / (Gamma(k + 1) * Gamma(n - 2 * k + 1));
            }

            return Gamma(n + 1) * x;
        }

        public static double LaguerreL(double n, double z)
        {
            if (n == 0) return (1.0);

            double L0 = 1.0;
            double L1 = 1.0 - z;
            for (int k = 1; k < n; k++)
            {
                double L2 = ((2 * k + 1 - z) * L1 - k * L0) / (k + 1);
                L0 = L1;
                L1 = L2;
            }
            return (L1);
        }

        public static double LaguerreL(double n, double a, double z)
        {
            double l1 = 0;
            double r = 1;

            for (int i = 1; i <= n; i++)
            {
                double y = l1;
                l1 = r;
                r = ((2 * i - 1 + a - z) * l1 - (i - 1 + a) * y) / i;
            }

            return r;
        }

        public static double LegendreP(double n, double z)
        {
            double x = 0;

            for (int k = 0; k <= Math.Floor(n / 2); k++)
            {
                x += Pow(-1, k) * Binomial(n, k) * Binomial(2 * n - 2 * k, n) * Pow(z, n - 2 * k);
            }

            return (1 / Pow(2, n)) * x;
        }

        public static double ChebyshevT(double n, double z)
        {
            if (n == 0) return (1.0);

            double T0 = 1.0;
            double T1 = z;
            for (int k = 1; k < n; k++)
            {
                double T2 = 2 * z * T1 - T0;
                T0 = T1;
                T1 = T2;
            }
            return (T1);
        }

        public static double ChebyshevU(double n, double z)
        {
            if(n == 0)
            {
                return 1;
            }
            if(n == 1)
            {
                return z;
            }

            double x = 0;

            for (int k = 0; k <= Math.Floor(n / 2); k++)
            {
                x += (Pow(-1, k) * Gamma(n - k + 1) * Pow(2 * z, n - 2 * k)) / (Fact(k) * Gamma(n - 2 * k + 1));
            }

            return x;
        }

        public static double GegenbauerC(double n, double a, double z)
        {
            double C0 = 1.0;
            if (n == 0) return (C0);

            double C1 = -2.0 * a * z;

            for (int k = 1; k < n; k++)
            {
                double C2 = (2.0 * z * (k + a) * C1 - (k - 1 + a) * C0) / (k + 1);
                C0 = C1;
                C1 = C2;
            }
            return (C1);
        }
        
        public static double JacobiP(double n, double a, double b, double z)
        {
            double x = 0;

            for (int k = 0; k <= n; k++)
            {
                 x += (Pochhammer(-n, k) * Pochhammer(a + b + n + 1, k) * Pochhammer(a + k + 1, n - k)) / Fact(k);
            }

            return (1 / Gamma(n + 1)) * x;
            //return 0;
        }

        public static double LegendreP(double n, double u, double cof, double z)
        {
            double x = 0;

            for (int k = 0; k <= n; k++)
            {
                x += ((Pochhammer(-n, k) * Pochhammer(n + 1, k)) / (Gamma(k - u + 1) * Fact(k))) * Pow((1 - z) / 2, k);
            }

            return (Pow(1 + z, u / 2) / Pow(1 - z, u / 2)) * x;
        }

        public static double LegendreP(double n, double u, double z)
        {
            return LegendreP(n, u, 2, z);
        }

        #endregion

        #region Fonctions harmoniques
        public static double Harmonic(double t)
        {
            if (Math.Truncate(t) != t) return 0;
            return Constant.EulerGamma + DiGamma(t + 1.0);
        }
        public static double GeneralHarmonic(double n, double m)
        {
            if (Math.Truncate(n) != n) return 0;
            double sum = 0;
            for (int i = 0; i < n; i++)
            {
                sum += Math.Pow(i + 1, -m);
            }

            return sum;
        }
        #endregion

        #region Autres fonctions utiles
        public static double Binomial(double n, double k)
        {
            return Gamma(n + 1) / (Gamma(k + 1) * Gamma(n - k + 1));
        }
        public static double Multinomial(params double[] a)
        {
            if(a.Any(x1 => Math.Truncate(x1) != x1)) throw new ParseException("", 0, "Tous les paramètres doivent être entiers");
            if(a.Any(x2 => x2 <= 0)) throw new ParseException("", 0, "Tous les paramètres doivent être strictement positifs.");

            double x = 1;
            a.All(y =>
            {
                x = x * Gamma(y + 1);
                return true;
            });

            return Gamma(a.Sum() + 1) / x;
        }
        public static double Pochhammer(double a, double n)
        {
            //if(Math.Truncate(-a) != -a) throw new ParseException((-a).ToString(), 0, "a doit être entier.");
            //if (!(-a >= 0)) throw new ParseException(a.ToString(), 0, "a doit être inférieur ou égal à 0");
            //if(!(a <= 0)) throw new ParseException(a.ToString(), 0, "a doit être inférieur ou égal à 0");
            //if (Math.Truncate(n) != n) throw new ParseException(n.ToString(), 0, "n doit être entier.");
            //if (!(n <= -a)) throw new ParseException(n.ToString(), 0, "n doit être inférieur ou égal à -a");


            //return Gamma(a + n) / Gamma(a);
            //return (Pow(-1, n) * Fact(-a)) / Fact(-a-n);
            //return (Pow(-1, n) * Gamma(-a + 1)) / Gamma(-a - n + 1);
            //return Gamma(a + 1) / Gamma(a - n + 1);
            return Binomial(a, n) * Gamma(n + 1);
        }
        #endregion

        #region Fonctions hypergéométriques
        public static double Hypergeometric1F1(double a, double b, double z)
        {
            double x = 0;
            double k = 1;

            while (true)
            {
                double nw = (Pochhammer(a, k) * Pow(z, k)) / (Pochhammer(b, k) * Fact(k));
                if (Math.Round(x + nw, 15) == Math.Round(x, 15))
                {
                    break;
                }
                x += nw;

                k++;
            }

            return x;
        }

        public static double ParabolicCylinderD(double v, double z)
        {
            return Pow(2, v / 2) * Math.Sqrt(Constant.Pi) * Pow(Constant.E, -(Pow(z, 2) / 4)) *
                   ((1 / Gamma((1 - v) / 2)) * Hypergeometric1F1(-v / 2, 1 / 2, z * z / 2) -
                    (Math.Sqrt(2) * z / Gamma(-v / 2)) * Hypergeometric1F1((1 - v) / 2, 3 / 2, z * z / 2));
        }
        #endregion

        #region Fonctions Bessel
        public static double BesselJ0(double x)
        {
            double ax;

            if ((ax = Math.Abs(x)) < 8.0)
            {
                double y = x * x;
                double ans1 = 57568490574.0 + y * (-13362590354.0 + y * (651619640.7
                    + y * (-11214424.18 + y * (77392.33017 + y * (-184.9052456)))));
                double ans2 = 57568490411.0 + y * (1029532985.0 + y * (9494680.718
                    + y * (59272.64853 + y * (267.8532712 + y * 1.0))));

                return ans1 / ans2;
            }
            else
            {
                double z = 8.0 / ax;
                double y = z * z;
                double xx = ax - 0.785398164;
                double ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4
                    + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
                double ans2 = -0.1562499995e-1 + y * (0.1430488765e-3
                    + y * (-0.6911147651e-5 + y * (0.7621095161e-6
                    - y * 0.934935152e-7)));

                return Math.Sqrt(0.636619772 / ax) *
                    (Math.Cos(xx) * ans1 - z * Math.Sin(xx) * ans2);
            }
        }
        public static double BesselJ(double x)
        {
            double ax;
            double y;
            double ans1, ans2;

            if ((ax = Math.Abs(x)) < 8.0)
            {
                y = x * x;
                ans1 = x * (72362614232.0 + y * (-7895059235.0 + y * (242396853.1
                    + y * (-2972611.439 + y * (15704.48260 + y * (-30.16036606))))));
                ans2 = 144725228442.0 + y * (2300535178.0 + y * (18583304.74
                    + y * (99447.43394 + y * (376.9991397 + y * 1.0))));
                return ans1 / ans2;
            }
            else
            {
                double z = 8.0 / ax;
                double xx = ax - 2.356194491;
                y = z * z;

                ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4
                    + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
                ans2 = 0.04687499995 + y * (-0.2002690873e-3
                    + y * (0.8449199096e-5 + y * (-0.88228987e-6
                    + y * 0.105787412e-6)));
                double ans = Math.Sqrt(0.636619772 / ax) *
                    (Math.Cos(xx) * ans1 - z * Math.Sin(xx) * ans2);
                if (x < 0.0) ans = -ans;
                return ans;
            }
        }

        public static double BesselJ(double nu, double x)
        {
            if(Math.Truncate(nu) != nu) throw new ParseException("", 0, "nu doit être entier.");
            int n = (int)Math.Truncate(nu);
            int j, m;
            double ax, bj, bjm, bjp, sum, tox, ans;
            bool jsum;

            double ACC = 40.0;
            double BIGNO = 1.0e+10;
            double BIGNI = 1.0e-10;

            if (n == 0) return BesselJ0(x);
            if (n == 1) return BesselJ(x);

            ax = Math.Abs(x);
            if (ax == 0.0) return 0.0;
            else if (ax > (double)n)
            {
                tox = 2.0 / ax;
                bjm = BesselJ0(ax);
                bj = BesselJ(ax);
                for (j = 1; j < n; j++)
                {
                    bjp = j * tox * bj - bjm;
                    bjm = bj;
                    bj = bjp;
                }
                ans = bj;
            }
            else
            {
                tox = 2.0 / ax;
                m = 2 * ((n + (int)Math.Sqrt(ACC * n)) / 2);
                jsum = false;
                bjp = ans = sum = 0.0;
                bj = 1.0;
                for (j = m; j > 0; j--)
                {
                    bjm = j * tox * bj - bjp;
                    bjp = bj;
                    bj = bjm;
                    if (Math.Abs(bj) > BIGNO)
                    {
                        bj *= BIGNI;
                        bjp *= BIGNI;
                        ans *= BIGNI;
                        sum *= BIGNI;
                    }
                    if (jsum) sum += bj;
                    jsum = !jsum;
                    if (j == n) ans = bjp;
                }
                sum = 2.0 * sum - bj;
                ans /= sum;
            }

            return x < 0.0 && n % 2 == 1 ? -ans : ans;
        }
        public static double BesselY0(double x)
        {
            if (x < 8.0)
            {
                double y = x * x;

                double ans1 = -2957821389.0 + y * (7062834065.0 + y * (-512359803.6
                    + y * (10879881.29 + y * (-86327.92757 + y * 228.4622733))));
                double ans2 = 40076544269.0 + y * (745249964.8 + y * (7189466.438
                    + y * (47447.26470 + y * (226.1030244 + y * 1.0))));

                return (ans1 / ans2) + 0.636619772 * BesselJ0(x) * Ln(x);
            }
            else
            {
                double z = 8.0 / x;
                double y = z * z;
                double xx = x - 0.785398164;

                double ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4
                    + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
                double ans2 = -0.1562499995e-1 + y * (0.1430488765e-3
                    + y * (-0.6911147651e-5 + y * (0.7621095161e-6
                    + y * (-0.934945152e-7))));
                return System.Math.Sqrt(0.636619772 / x) *
                    (System.Math.Sin(xx) * ans1 + z * System.Math.Cos(xx) * ans2);
            }
        }
        public static double BesselY(double x)
        {
            if (x < 8.0)
            {
                double y = x * x;
                double ans1 = x * (-0.4900604943e13 + y * (0.1275274390e13
                    + y * (-0.5153438139e11 + y * (0.7349264551e9
                    + y * (-0.4237922726e7 + y * 0.8511937935e4)))));
                double ans2 = 0.2499580570e14 + y * (0.4244419664e12
                    + y * (0.3733650367e10 + y * (0.2245904002e8
                    + y * (0.1020426050e6 + y * (0.3549632885e3 + y)))));
                return (ans1 / ans2) + 0.636619772 * (BesselJ(x) * Ln(x) - 1.0 / x);
            }
            else
            {
                double z = 8.0 / x;
                double y = z * z;
                double xx = x - 2.356194491;
                double ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4
                    + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
                double ans2 = 0.04687499995 + y * (-0.2002690873e-3
                    + y * (0.8449199096e-5 + y * (-0.88228987e-6
                    + y * 0.105787412e-6)));
                return System.Math.Sqrt(0.636619772 / x) *
                    (System.Math.Sin(xx) * ans1 + z * System.Math.Cos(xx) * ans2);
            }
        }
        public static double BesselY(double nu, double x)
        {
            if (Math.Truncate(nu) != nu) throw new ParseException("", 0, "nu doit être entier.");
            int n = (int)Math.Truncate(nu);
            double by, bym, byp, tox;

            if (n == 0) return BesselY0(x);
            if (n == 1) return BesselY(x);

            tox = 2.0 / x;
            by = BesselY(x);
            bym = BesselY0(x);
            for (int j = 1; j < n; j++)
            {
                byp = j * tox * by - bym;
                bym = by;
                by = byp;
            }
            return by;
        }
        public static double BesselI0(double x)
        {
            double ans;
            double ax = Math.Abs(x);

            if (ax < 3.75)
            {
                double y = x / 3.75;
                y = y * y;
                ans = 1.0 + y * (3.5156229 + y * (3.0899424 + y * (1.2067492
                   + y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
            }
            else
            {
                double y = 3.75 / ax;
                ans = (Math.Exp(ax) / Math.Sqrt(ax)) * (0.39894228 + y * (0.1328592e-1
                   + y * (0.225319e-2 + y * (-0.157565e-2 + y * (0.916281e-2
                   + y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1
                   + y * 0.392377e-2))))))));
            }

            return ans;
        }
        public static double BesselI(double x)
        {
            double ans;

            double ax = Math.Abs(x);

            if (ax < 3.75)
            {
                double y = x / 3.75;
                y = y * y;
                ans = ax * (0.5 + y * (0.87890594 + y * (0.51498869 + y * (0.15084934
                   + y * (0.2658733e-1 + y * (0.301532e-2 + y * 0.32411e-3))))));
            }
            else
            {
                double y = 3.75 / ax;
                ans = 0.2282967e-1 + y * (-0.2895312e-1 + y * (0.1787654e-1 - y * 0.420059e-2));
                ans = 0.39894228 + y * (-0.3988024e-1 + y * (-0.362018e-2 + y * (0.163801e-2 + y * (-0.1031555e-1 + y * ans))));
                ans *= Math.Exp(ax) / Math.Sqrt(ax);
            }
            return x < 0.0 ? -ans : ans;
        }
        public static double BesselI(double nu, double x)
        {
            if (Math.Truncate(nu) != nu) throw new ParseException("", 0, "nu doit être entier.");
            int n = (int)Math.Truncate(nu);

            if (n < 0)
                throw new ParseException("", 0, "n doit être positif.");
            else if (n == 0)
                return BesselI0(x);
            else if (n == 1)
                return BesselI(x);

            if (x == 0.0)
                return 0.0;

            double ACC = 40.0;
            double BIGNO = 1.0e+10;
            double BIGNI = 1.0e-10;

            double tox = 2.0 / Math.Abs(x);
            double bip = 0, ans = 0.0;
            double bi = 1.0;

            for (int j = 2 * (n + (int)Math.Sqrt(ACC * n)); j > 0; j--)
            {
                double bim = bip + j * tox * bi;
                bip = bi;
                bi = bim;

                if (Math.Abs(bi) > BIGNO)
                {
                    ans *= BIGNI;
                    bi *= BIGNI;
                    bip *= BIGNI;
                }

                if (j == n)
                    ans = bip;
            }

            ans *= BesselI0(x) / bi;
            return x < 0.0 && n % 2 == 1 ? -ans : ans;
        }
        #endregion

        #region Fonctions Zeta
        public static double Zeta(double s)
        {
            double x = 0;
            double k = 1;

            while (true)
            {
                double nw = 1 / Pow(k, s);
                if (Math.Round(x + nw, 15) == Math.Round(x, 15))
                {
                    break;
                }
                x += nw;

                k++;
            }

            return x;
        }
        public static double Zeta(double s, double a)
        {
            double x = 0;
            double k = 0;

            while (true)
            {
                double nw = 1 / Pow(Pow(a + k, 2), s / 2);
                if (Math.Round(x + nw, 15) == Math.Round(x, 15))
                {
                    break;
                }
                x += nw;

                k++;
            }

            return x;
        }
        public static double LerchPhi(double z, double s, double a)
        {
            double x = 0;
            double k = 0;

            while (true)
            {
                double nw = Pow(z, k) / Pow(Pow(a + k, 2), s / 2);
                if (Math.Round(x + nw, 15) == Math.Round(x, 15))
                {
                    break;
                }
                x += nw;

                k++;
            }

            return x;
        }
        #endregion

        #region Fonctions polylogarithmiques
        public static double PolyLog(double v, double z)
        {
            double x = 0;
            double k = 1;

            while (true)
            {
                double nw = Pow(z, k) / Pow(k, v);
                if (Math.Round(x + nw, 15) == Math.Round(x, 15))
                {
                    break;
                }
                x += nw;

                k++;
            }

            return x;
        }
        
        #endregion

        static Function()
	    {
	        FunctionsInf = new Dictionary<string, FunctionCallbackInf>();
            Functions4 = new Dictionary<string, FunctionCallback4>();
	        Functions3 = new Dictionary<string, FunctionCallback3>();
	        Functions2 = new Dictionary<string, FunctionCallback2>();
	        Functions = new Dictionary<string, FunctionCallback>();
	        FunctionsNone = new Dictionary<string, FunctionCallbackNone>();

	        Functions["abs"] = Math.Abs;
	        Functions["fact"] = Fact;
	        Functions["fact2"] = Fact2;

	        Functions["sin"] = Sin;
	        Functions["arcsin"] = Math.Asin;
	        Functions["sinh"] = Sinh;
	        Functions["arcsinh"] = Arcsinh;
	        Functions["csc"] = Csc;
	        Functions["arccsc"] = Arccsc;
	        Functions["csch"] = Csch;
	        Functions["arccsch"] = Arccsch;

	        Functions["sinc"] = Sinc;

	        Functions["cos"] = Cos;
	        Functions["arccos"] = Math.Acos;
	        Functions["cosh"] = Cosh;
	        Functions["arccosh"] = Arccosh;
	        Functions["sec"] = Sec;
	        Functions["arcsec"] = Arcsec;
	        Functions["sech"] = Sech;
	        Functions["arcsech"] = Arcsech;

	        Functions["tan"] = Tan;
	        Functions["arctan"] = Math.Atan;
	        Functions2["arctan"] = Math.Atan2;
	        Functions["tanh"] = Tanh;
	        Functions["arctanh"] = Arctanh;
	        Functions["cot"] = Cot;
	        Functions["arccot"] = Arccot;
	        Functions["coth"] = Coth;
	        Functions["arccoth"] = Arccoth;

	        Functions["sqrt"] = Math.Sqrt;
	        Functions["curt"] = CubeRoot;
	        Functions2["root"] = NthRoot;
	        Functions["log10"] = Math.Log10;
	        Functions2["log"] = Math.Log;
	        Functions["ln"] = Ln;
	        Functions["exp"] = Math.Exp;

	        Functions["ceil"] = Math.Ceiling;
	        Functions["floor"] = Math.Floor;
	        FunctionsInf["gcd"] = GCD;
	        FunctionsInf["lcm"] = LCM;
	        FunctionsInf["max"] = Max;
	        FunctionsInf["min"] = Min;

	        FunctionsInf["polynomial"] = Polynomial;
	        Functions2["binomial"] = Binomial;
	        FunctionsInf["multinomial"] = Multinomial;
	        Functions2["hermiteh"] = HermiteH;
	        Functions2["laguerrel"] = LaguerreL;
            Functions3["laguerrel"] = LaguerreL;
	        Functions2["legendrep"] = LegendreP;
	        Functions2["chebyshevt"] = ChebyshevT;
	        Functions2["chebyshevu"] = ChebyshevU;
	        Functions3["gegenbauerc"] = GegenbauerC;
	        Functions4["jacobip"] = JacobiP;
	        Functions3["legendrep"] = LegendreP;
	        Functions4["legendrep"] = LegendreP;

            // fonctions bessel
            Functions["besselj0"] = BesselJ0;
            Functions["besselj"] = BesselJ;
            Functions2["besselj"] = BesselJ;
            Functions["bessely0"] = BesselY0;
            Functions["bessely"] = BesselY;
            Functions2["bessely"] = BesselY;
            Functions["besseli0"] = BesselI0;
            Functions["besseli"] = BesselI;
            Functions2["besseli"] = BesselI;

	        FunctionsInf["kroneckerdelta"] = KroneckerDelta;
	        Functions["kroneckerdelta"] = KroneckerDelta;

	        // Fonctions Beta
	        Functions2["betaln"] = BetaLn;
	        Functions2["beta"] = Beta;
	        Functions3["betainc"] = BetaIncomplete;
	        Functions3["betareg"] = BetaRegularized;


	        // Fonctions Gamma
	        Functions["gamma"] = Gamma;
	        Functions["gammaln"] = GammaLn;
	        Functions2["gammaupreg"] = GammaUpperRegularized;
	        Functions2["gammaupinc"] = GammaUpperIncomplete;
	        Functions2["gammalowinc"] = GammaLowerIncomplete;
	        Functions2["gammalowreg"] = GammaLowerRegularized;
	        Functions2["gammalowreginv"] = GammaLowerRegularizedInv;
	        Functions["digamma"] = DiGamma;
            Functions["polygamma"] = PolyGamma;
	        Functions["digammainv"] = DiGammaInv;
	        Functions["loggamma"] = LogGamma;

            // Fonctions Zeta
            Functions["zeta"] = Zeta;
            Functions2["zeta"] = Zeta;
            Functions3["lerchphi"] = LerchPhi;

            // Fonctions polylogarithmiques
            Functions2["polylog"] = PolyLog;

	        // Fonctions Erf
	        Functions["erf"] = Erf;
	        Functions["erfc"] = Erfc;
	        Functions["erfinv"] = ErfInv;
	        Functions2["erfimpl"] = ErfImp;
	        Functions3["erfimplinv"] = ErfInvImpl;

            Functions2["paraboliccylinderd"] = ParabolicCylinderD;
            Functions3["hypergeometric1f1"] = Hypergeometric1F1;

	        Functions["harmonic"] = Harmonic;
	        Functions2["genharmonic"] = GeneralHarmonic;


	        Functions["int"] = Math.Truncate;
	        Functions["fibo"] = Fibo;
	        Functions2["rnd"] = Round;
	        Functions2["eq"] = Eq;
            Functions["lucasl"] = LucasL;
            Functions["productlog"] = ProductLog;
            Functions2["productlog"] = ProductLog;

	        FunctionsInf["avr"] = Average;
            FunctionsInf["med"] = Mediane;
            FunctionsInf["q1"] = Q1;
            FunctionsInf["q3"] = Q3;
            FunctionsInf["sum"] = Somme;
	        FunctionsNone["rand"] = Rand;
	        Functions2["rand"] = Rand;
	        Functions2["randn"] = RandN;


	        FunctionsList = new List<string>
	        {
	            "pi",
	            "e",
	            "phi",
	            "eulergamma",
	            "catalan",
	            "c"
	        };
            Functions.All(x =>
            {
                FunctionsList.Add(x.Key + "(a)");
                return true;
            });
            Functions2.All(x =>
            {
                FunctionsList.Add(x.Key + "(a; b)");
                return true;
            });
            Functions3.All(x =>
            {
                FunctionsList.Add(x.Key + "(a; b; c)");
                return true;
            });
            Functions4.All(x =>
            {
                FunctionsList.Add(x.Key + "(a; b; c; d)");
                return true;
            });
            FunctionsInf.All(x =>
            {
                FunctionsList.Add(x.Key + (x.Key == "polynomial" ? "(z; a; b; c; ...)" : "(a; b; c; ...)"));
                return true;
            });
	    }
		public override bool Signed
		{
			get { return m_signed; }
			set { m_signed = value; }
		}

		public override double Value
		{
			get
			{
			    if (FunctionsInf.ContainsKey(m_func) && m_terms.Count > 1)
			    {
                    var value = FunctionsInf[m_func](m_terms.Select(x => x.Value).ToArray());
                    if (Signed)
                        return -value;
                    return value;
			    }
			    else
			    {
			        if (m_terms[0].Stack.Count == 0)
			        {
			            var value = FunctionsNone[m_func]();
			            if (Signed)
			                return -value;
			            return value;
			        }
			        else if (m_terms.Count == 1)
			        {
			            var value = Functions[m_func](m_terms[0].Value);
			            if (Signed)
			                return -value;
			            return value;
			        }
			        else if (m_terms.Count == 2)
			        {
			            var value = Functions2[m_func](m_terms[0].Value, m_terms[1].Value);
			            if (Signed)
			                return -value;
			            return value;
			        }
			        else if (m_terms.Count == 3)
			        {
			            var value = Functions3[m_func](m_terms[0].Value, m_terms[1].Value, m_terms[2].Value);
			            if (Signed)
			                return -value;
			            return value;
			        }
                    else if (m_terms.Count == 4)
                    {
                        var value = Functions4[m_func](m_terms[0].Value, m_terms[1].Value, m_terms[2].Value, m_terms[3].Value);
                        if (Signed)
                            return -value;
                        return value;
                    }
			    }

                throw new ParseException(string.Empty, 0, string.Format("Aucune fonction acceptant {0} arguments", m_terms.Count));
			}
		}
		static public bool IsValidDigit(string equation, int index)
		{
			var spi = Term.ExtractName(equation, index).ToLower();
			if (spi.Length < 2)
				return false;
			if (Functions.ContainsKey(spi))
				return true;
			if (Functions2.ContainsKey(spi))
				return true;
            if (Functions3.ContainsKey(spi))
                return true;
		    if (Functions4.ContainsKey(spi))
		        return true;
            if (FunctionsInf.ContainsKey(spi))
                return true;
			if (FunctionsNone.ContainsKey(spi))
				return true;
			return false;
		}
		public override bool IsValid(string equation, int index)
		{
			return IsValidDigit(equation, index);
		}
		
		string m_func = string.Empty;
		int m_startIndex;
		int m_endindex;
		List<Term> m_terms;
		public override int Length
		{
			get
			{
				return m_endindex- m_startIndex + 1;
			}
		}
		public Function(string src, int start, EquationElement root)
		{
			m_startIndex = start;
			m_func = Term.ExtractName(src, start).ToLower();
			start += m_func.Length;
			// sauter tous les espaces, mais le prochain caractère doit être une parenthèse ouvrante
			while (src[start] == ' ')
				start++;
			if (src[start] != '(')
				throw new ParseException(src, m_startIndex, "La fonction doit commencer par '('");
			var termstart = start;
			var end = Term.FindMatchingEnd('(', src, termstart);
			if (end < termstart)
				throw new ParseException(src, m_startIndex, "Pas de parenthèse fermante correspondante trouvée");

			m_endindex = end;
			var allterms = src.Substring(termstart+1, end-termstart-1);
			//string[] terms = allterms.Split(',');
			var terms = GetTerms(allterms);
			m_terms = new List<Term>();
			foreach (var term in terms)
			{
				var newterm = new Term();
				newterm.Parse(term, root);
				newterm.Parent = this;
				m_terms.Add(newterm);
			}
		}

		private IEnumerable<string> GetTerms(string allterms) 
		{ 
			var splitIndex = allterms.IndexOf(';');
			if (splitIndex < 0)
				return new[] {allterms};
			
			var startChars = new[] { '(','[', ';' };
			var start = 0;
			var result = new List<string>();
			while (start < allterms.Length)
			{
				var ss = allterms.IndexOfAny(startChars, start);
				if (ss < 0)
				{
					var s = allterms.Substring(start).Trim();
					if (s.Length > 0)
						result.Add(s);
					break;
				}
				int len;
				if (allterms[ss] == ';')
				{
					// copier depuis le début
					len = ss - start;
					var s = allterms.Substring(start, len).Trim();
					if (s.Length > 0)
						result.Add(s);
					start = ss+1;
					continue;
				}
				var termend = Term.FindMatchingEnd(allterms[ss], allterms, ss);
				len = termend - start+1;
				var su = allterms.Substring(start, len).Trim();
				if (su.Length > 0)
					result.Add(su);
				start = termend+1;
				ss = allterms.IndexOfAny(startChars, start);
				if (ss > 0)
					start = ss+1;
			}
			return result.ToArray();
		}


		public override string ToString()
		{
			var sb = new StringBuilder();
			sb.Append(m_func);
			sb.Append('(');
			for (var index = 0; index < m_terms.Count; index++)
			{
				sb.Append(m_terms[index]);
				if (index < m_terms.Count-1)
					sb.Append(';');
			}
			sb.Append(')');
			return sb.ToString();
		}
	}
}
